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We discuss t-J-U model on a honeycomb monolayer that has the same low-energy description of 
the kinetic term as graphene bilayer, and in particular study coexistence of antiferromagnetism and 
superconducting correlations that originate from Cooper pairs without phase coherence. We show 
that the model is relevant for the description of graphene bilayer and that the presence of the d + id 
superconducting correlations with antiferromagnetism can lead to quadratic dependence in small 
magnetic fields of the gap of the effective monolayer consistent with the transport measurements of 
Velasco et al. on the graphene bilayer. 
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I. INTRODUCTION 

The interaction effects are important for the physics 
of graphene bilayer; recent experiments reveal gapped 
phase(s) in the undoped graphene bilayer which without 
interactions would represent a gapless system. In a re- 
cent experiment, in Ref. [ll, with high quality samples, 
a completely insulating behavior was detected in trans- 
port measurements. Theoretical investigations, mean 
field and renormalization group approaches [2hl6j , speak 
for a close competition of a few, mostly gapped, phases. 
One of the most prominent candidates for an explana- 
tion of the experiment in Ref. [l] is a layer antiferromag- 
netic (LAF) state. The main reason for the existence of 
this state would be an on-site Coulomb repulsion, U; in- 
deed as pointed out in Ref. ,6i, a Hubbard model on a 
honeycomb bilayer lattice would lead to the LAF state, 
both in weak and strong U limit. This may remind us of 
the behavior of the Hubbard model on the square lattice 
and the antiferromagnetic (AF) behavior due to nesting 
in the weak coupling limit. The Hubbard model on a 
square lattice is usually invoked as a model for cuprates 
in its strong coupling limit which forbids the double oc- 
cupancy and leads to a "perfect" AF behavior at half 
filling. On the other hand the estimate for U is hard 
to know in the graphene bilayer and certainly depends 
on the computational scheme but it is expected to be 
stronger than both (inter and intra-layer) hoppings. Due 
to the smallness of the gap revealed in the experiment in 
Ref. [l| we will not consider the large U limit (exclusion 
of double occupancy) when modeling graphene bilayer. 
But we will keep the on-site repulsion as a main cause of 
the insulating behavior detected in the experiment. As 
expected from previous approaches this will lead to AF 
insulating behavior but seems not sufficient to describe 
all phenomena detected in the experiment. An additional 
order parameter, besides the one that describes the anti- 
ferromagnetism is necessary for the complete explanation 
of the transport data of the experiment [T], Q . 

In this work we will look for the additional order pa- 



rameter that can coexist with antiferromagnetism in the 
graphene bilayer at half filling. We will argue that this is 
d + id (broken time reversal symmetry) - wave supercon- 
ducting order parameter. This {d + id) order parameter 
and its coexistence with antiferromagnetism was already 
found at finite (non-zero) dopings in a numerical (Grass- 
man tensor product) approach to t — J (large U) model 
on the honeycomb monolayer in Ref. jTJ Due to the as- 
sumed moderate (not large) value of U in our model of 
graphene bilayer the AF and d + id superconducting or- 
der parameter can coexist even at half-filling. Our model 
of the graphene bilayer can be described as a t-J-U model 
on an effective honeycomb lattice and, in the following, 
we will argue why this model is relevant for the descrip- 
tion of graphene bilayer. 

II. MODEL AND ITS MOTIVATION 

The kinetic part of the Hamiltonian that describes the 
garphene bilayer on two honeycomb lattices, which are 
Bernal stacked, is 

Hn = —t y^(a! - 6, - , ? + al ^ b„ ^ ; + h.c.) 

/ ^ / ^ ^ l,n,cr l,n+d,a 2.n,a 2,n — d.cr / 
n,a s 

+UXl("lfl,'T"2,H,<T + /l-C). (1) 
n,(T 

The index i = 1,2 denotes the layer index. In Fig. 1 the 
relative positions of two triangular sublattices, Ai and 
Bi, for the lattice 1, and A2 and B2, for the lattice 2 are 
illustrated. In Eq.(IT]) t is the hopping energy between 
nearest sites in each layer, and t± is the same energy for 
hopping between the layers. The on-site creation (anni- 
hilation) operators, a]^ - „{a,i,n,a), are for the electrons in 
the sublattice Ai of the layer i with spin a =t, i, and 
6| - ^(6i.fj_cr) for the electrons in the sublattice Bi. S''s 
are defined as 61 — a(0, l/-\/3), S2 — a/2(l, —1/^/3), and 
S3 — a/2(— 1, — 1/-\/3), and a — VS Ucc, dec is the dis- 
tance between sites and a is the next to nearest neighbor 
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FIG. 1: A view of Bernal stacked honeycomb lattices 1 and 
2 with corresponding sublattice sites Al and Bl, and A2 and 
B2. 



distance. 

Accordin g; to the results of renormalization group ap- 
proaches [5l4lOl| the antiferromagnetic (LAF) state and a 
nematic state are main instabilities that may arise due to 
interactions in the graphene bilayer. The nematic state 
is a result of an ordering in the particle-hole channel that 
can be described by a non-zero expectation value of the 
hopping (bond) operator, 



/-^ l,ri+i5i, 



r 2,n — do 



(2) 



where n + Si and n— 82 denote sites that are near neigh- 
bors on the honeycomb lattice that make sublattices Bl 
and B2, as shown in Fig. 1. On the same lattice the an- 
tiferromagnetic ordering can occur which describes the 
LAF state. 

The nematic i.e. bond ordering as described in Ref. 7 
can be thought as a dx^-y2 CDW (or a d + id CDW as 
described in Ref. |^ - see below Eq. for the definition 
of the d^2_y2 and d + id ordering) in the language of 
Ref. O In the same reference this instability ("hidden 
order" on a square lattice) was proposed for cuprates. 
Just as in the case of the square lattice and cuprates, 
this ordering, in the graphene bilayer, can be thought as 
the result of short-range interactions and superexchange 
processes. These interactions and processes can also lead 
to the LAF state. The effective Heisenberg interaction 
on the Bl and B2 honeycomb, of neighboring sites i and 
j (ri + S and n ~ 6') can be rewritten in terms of the 
hopping operators (Eq.©) as 



SiSj 



2ja 



1 

4' 



Thus, on phenomenological grounds, we will assume that 
the value of J is an independent parameter in the effec- 
tive description, favorable for both instabilities, nematic 
and antiferromagnetic. 

With respect to cuprates the on-site repulsion, U, is 
not that strong to preclude the double occupancy on the 
graphene bilayer. Therefore we will explicitly include this 
interaction in our model Hamiltonian, which can be de- 
scribed as 



(4) 



where 



(3) 



-./^ E ^fci(n + ^i) 5'fc2(r^-^2), (5) 

" .51,52 

and the summation is over the near neighbors, and 

kelBZ i=l,2 Hi 

with 7fc = ^^exp{ifc(5}, fii = fi + Si and n2 = n — 62- 
The kinetic term (written above in the momentum space) 
in real space, on the effective honeycomb monolayer 
of Bl and B2 sites, describes near-neighbor and two 
times weaker third-neighbor hopping. It can be recov- 
ered in the case of the noninteracting honeycomb bilayer 

by taking interlayer hopping to be large. In that case 
.2 

tcs = j^- In the small momentum limit the kinetic term 
in Eq.(l6]) becomes the one of the graphene bilayer, i.e. 
teS ll j^{kx T ikyf near K points: K± = ±^(|,0) 
[19I [20I I . Because we look for the low energy properties 
we will keep the two-band extension with 7^ throughout 
the Brillouin zone. 

This model is similar to the t — J — U model defined 
on the square lattice and known from previous investiga- 
tions. The t — J — U model appeared also in the context 
of gossamer superconductivity [21.-25] , the superconduc- 
tivity that can exist even at half-filling. 

III. MEAN FIELD APPROACH 

In order to apply a mean field approach we can use the 
identity, 

S^S, = -lib^b,, h,b,,){blbl, - b]^bl) + i. (7) 
We will define the superconducting order parameter, 

A5- = {butb^.+si - ^2.+5t^i4) > (8) 



and should be a part of an effective description of the 
graphene bilayer. According to Refs. andfiol the only 
instability in the weak coupling limit of the Hubbard 
model on the graphene bilayer is antiferromagnetism. 



where S can be any of the near-neighbor vectors on the 
honeycomb lattice, and 



m = {riif - nil), 



(9) 
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the antiferromagnetic order parameter. In the fohowing 
we win use the fohowing notation, t = teff j and apply to 
the Hamiltonian in Eq.Q the mean field ansatzes. We 
will generalize the derivation of Ref. '26' to the case with 
the antiferromagnetic order parameter. 
If we use spinors, 



we can write the mean field Hamiltonian as 



(10) 



n^, + lNjJ2\^f + U^N, (11) 



where N is the number of unit cells, 







-til 







n = 


-til' 


-|m 
-4A*(-fc) 


-^A(-fc) 










till 




_ -4A*(fc) 





tl-k 





(12) 

and A(fc) — A^expjiM}. The symmetry analysis of 
the order parameter on a honeycomb lattice, first done 
in Ref. [27|, concluded that there are three possibilities. 



A(l, 1, 1) 
A (2, -1, -1) 
A(0, 1, -1) 



(13) 



that span the space of order parameter. The last two 
possibilities belong to a two-dimensional subspace of ir- 
reducible representation of permutation group (28| . 
The s-wave, Ag = A(l, 1, 1), has nodes at K points be- 
cause Afc — A'jk- For dx2_y2 wave. As = A(2,— 1,— 1), 



= A[3± ^(fc, ± 
1), the expansions 



near K± points we have Ai{K± + k) 
iky)], and for dxy wave, A5 — A(0, 1, - 
are A2{K± + k) = Ai[±^/?, - ^{kx ± iky)]. Therefore if 
the d + id combination, Ai(fc) ± i\/3A2{k) is taken near 
one of the K points the order parameter is a constant 
(6 A) and at the other K point is linear in kx and ky. 
Therefore, instead of having the coefficients of the same 
absolute magnitude with &Jj,-f-62~fe4. ^^'^ ^2fct^i-fc4. (^^'i 
b2-kibikf and 6i-fc4,&2fet) for a fixed valley point d + id 
singles out one spin projection (up or down) to be asso- 
ciated with sites on layer 1 and the opposite one to be 
associated with sites on layer 2. Thus it favors pairing 
(Cooper pairs) in which layer index is associated with def- 
inite spin projection just as in an antiferromagnetically 
ordered state i.e. LAF state described above. 

d + id wave and s wave can coexist, as rotationally 
symmetric states, with the LAF state although only for 
certain values of J and U parameters. One can show that 
for the coexistence of the LAF state and s-wave J ^ U, 
which cannot be the case in the graphene bilayer. For 
d+id wave, on the other hand, one can find an interval for 
couplings, J and U, for which the LAF state and d + id- 
wave pairing can coexist. Expanding the mean field equa- 
tions to fourth order in the ratio of superconducting and 



antiferromagnetic order parameter and in the weak cou- 
pling limit i.e. t > J,U we find | — jJw < < | + jJw, 
We expect that the interval will broaden 



with w 



when the short-range correlations (due to U) are prop- 
erly taken into account that will renormalize (reduce) the 
effective value of t. This was worked out for square lat- 
tice in Ref. "23 in a comprehensive (renormalized mean 
field) study oi t — J — U model, and a interval of cou- 
plings was identified for which antiferromagnetism and 
superconducting correlations can coexist at half-filling. 
Furthermore, in Ref. 2.4 variational Monte Carlo method 
was applied to the same system, and a finite value of the 
pairing amplitude (A) was found in the antiferromag- 
netic region (with no superconducting phase coherence). 
We expect a similar situation in our case. 

As in the case of honeycomb monolayer in Ref. [27l we 
can diagonalize the free part of the above Hamiltonian 
and come to the following expressions for order parame- 
ters. 



and 



where 



arg{j{k)). Due to the expansion of A;; 



around K points in the case of d + id wave we have 



C. 



K± + k ^K±+k 



[kx i iky)'^ 



(14) 



where the last sign is independent of K points. Thus we 
recovered the basic signatures of d -I- id pairing; (a) the 
order parameter is an eigenfunction of orbital angular 
momentum with eigenvalue equal to two, and (b) due 
to the same sign (chirality) at both K points this wave 
is time reversal symmetry breaking wave on the bilayer 
honeycomb lattice. Therefore in analyzing this wave we 
can keep the leading behavior in Ai(A:) and A2(A:) as this 
effectively captures the basic phenomenology of d -f id- 
wave. Thus we take 





r fm 


-til 





-J3A" 




-til' 










id — 










till 




_ -J3A 





tl'-k 


-f m _ 



(15) 



in the case of A(fc) — Ai(fc) — i\/iA2{k) combination, or 
-til 



n 



d-\-id 



-ti 



*2 





J3A 



-J3A 
t-fl 






ti*J 



(16) 



in the case when A(fc) — Ai(fc) +i\/3A2(fc). In the same 
low-momentum limit 'yK±+k ~ To,^{kx T iky). We will 
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m, J3A = A, and ta^j = t in and V denotes the perturbation, 



the following. 

We take A to be purely real and without t he p hase 
(U(l)) degree of freedom i.e. phase coherence [23] that 
would lead to supercurrents proportional to the gradi- 
ent of this phase that would screen a magnetic field that 
may be present. We assume that supercurrents cannot 
develop in the antiferromagnetic, insulating background. 

The Bogoliubov spectrum is the same irrespective 
whether we ask for energy eigenvalues in the case defined 
by Eq. pS)) or Ea.([TC|. and with the introduced redefini- 
tions the eigenvalues are 



E = ±] 



,A 



A. 



+ (-± Jt2fc4 + (-)2) 



(17) 



Therefore the two different chirality states of the d-wave 
are equally likely in the presence of the antiferromagnetic 
ordering. 



IV. PRESENCE OF SMALL MAGNETIC FIELD 

In the presence of magnetic field, due to the minimal 
prescription, we may introduce a pair of creation and 
annihilation operators and express the resulting Hamil- 
tonian matrix around if^ point as 



HI 



{d-id)B 



m 



-A 



t^2 






TO 



-A 


c^c(at)2 



— TO 




(18) 

Here we introduced Wc = — , where B is the magnetic 
field and to is the effective mass of the graphene bilayer, 
^ = 2t. The eigenvectors can be expressed as 4-spinor 
coefficients of eigenvectors of a^a - operator, a^a'^n = 
n^„, classified by integer eigenvalues n: 0,1,2,.... In 
the presence of small magnetic field we will look for the 
eigenstates in the form 



*n = [Cl C2 C3 CaY 



0,1,2, 



(19) 



The Nambu-Gorkov formalism with 4-spinors artifi- 
cially doubles the degrees of freedom. This appears in 
spectra as doubling of energy levels {^E). Thus when 
solving the H+(jj_id)B we have to keep levels that are 
continuously related to energy levels with no supercon- 
ducting instability (A ^ 0) and are pertinent to the 2x2 
upper, left block of the Hamiltonian matrix. 

The Hamiltonian in Eq. (|T8| we will consider under the 
approximation of a small magnetic field and rewrite it as 



'^°+(d-id)B — ^0 



V 



where 



Hn 



m -A 

-m 

TO 

-A -m 



(20) 



(21) 



V = 





-ujc{a) 





-ujc[a 






t^2 
















(22) 



Taking as solutions only values that are connected con- 
tinuously in the limit A — >■ to the upper 2x2 left part 
of Hq we get for the eigenvalues and eigenvectors of Hq: 



E'^ 



"A^ 



[0, 1, 0, 0]^|n), 
c[m + E, 0, 0, ~A]' 



where E — \/r. 



A^ and c 



Considering 



y/2E{E+'n 

the small magnetic field to second order as perturbation 



we get El = — m- 



(n+2)(K+l) uj^ 



and EI^ = E+ 



n(n — 1) Lj^ 



Considering the same problem at K' = —K point we get 
E^ = m+ ("+2K"+i) ^ and E^ ^ -E - 

Thus, by analyzing the spectra of both K points to- 
gether, we can conclude that with the inclusion of small 
magnetic fields the gap changes from Eg — 2m value to 

2 

Eg — 2m -f 2-^ in the presence of d — id correlations. 
Without the correlations or with d + id correlations the 
gap will not have the correction quadratic in small mag- 
netic field, which direction is fixed in Ea.([T8|. d — id 
correlations minimize the energy of the system by shift- 
ing also the energy levels closest to the Fermi point. In 
the Appendix we compare the energies of the states with 
d + id and d — id correlations, and show that d — id are 
indeed of the lower energy. 

The energy minimization, when the direction of per- 
pendicular magnetic field is opposite requires that the 
superconducting correlations are of d 4- id kind. Thus the 
change in the direction of magnetic field is followed by the 
change in the chirality of the superconducting instability 
i.e. B — )■ —B followed by d — id — >■ d -f id as one would 
expect from the superconducting instability that has or- 
bital and therefore magnetic moment. This amounts to 
just switching of the previously found spectra between K 
and K' points. The gap is the same irrespective of the 
direction of the magnetic field although with the inclu- 
sion of superconducting correlations linear in k (on the 
diagonal in the Eq.lfTS]) and Eq. ([T6| ) leads to an asym- 
metry which may be related to the asymmetry detected 
in the transport measurements of Ref. X with respect to 
the change in the direction of the external field. 



V. DISCUSSION AND CONCLUSION 

In the literature we find several proposals, Refs. Isol - 
[32l . for the explanation of the data of Ref. iH See also 
Ref. [3^ for further experimental investigations on the 
same system and possible explanations based on anoma- 
lous quantum Hall physics. Ref. Isflj by Kharitonov in- 
troduces an additional order parameter to the Neel order 
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parameter, but the resulting gap dependence does not 
have a minimum at B (magnetic field) =0 - compare Fig. 
3 in Ref. HO, in contrast to what can be seen from the 
transport measurements - compare Fig. 3 in Ref. T". Ref. 
Isil . by Zhu, Aji, and Varma, with the interesting pro- 
posal of taking into account the full four band structure, 
still gives linear dependence on B of the gap - compare 
with Fig. 6 in Ref. [3l| in contrast to the quadratic de- 
pendence on small B as seen in the experiment. Ref. 
|32 by Roy does describe the quadratic dependence based 
on a mean field treatment of spin magnetism, where a 
phenomenological ferromagnetic interaction next to the 
Neel ordering among the spins of electrons is introduced 
which existence (with a precise magnitude) is necessary 
to obtain a correspondence to the experimental data. 

Our approach is also mean field and phenomenological, 
though clearly motivated microscopically by the physics 
of the t — J — U model, as we introduce orbital mag- 
netism of superconducting correlations. This leads to the 
quadratic dependence of the gap on small B as observed 
in the experiment. Thus we demonstrated a possibility 
that the quadratic dependence on small magnetic field 
observed in the experiment of Ref. [l| may be due to the 
time reversal symmetry breaking d-wave superconduct- 
ing correlations that coexist with antiferromagnetism. 

The d + id wave superconductivity and antiferromag- 
netism at high dopings of the graphene monolayer was 
studied in Refs. Q and IH. It was shown [IH that both 
instabilities are connected with the existence of the on- 
site repulsion, U. 

The last and important question we would like to dis- 
cuss is how our proposal can explain the behavior of the 
system gap in strong (and moderate) magnetic fields. In 
other words the question is how does the antiferromag- 
netic ground state with d + id superconducting correla- 
tions evolve in the many-body state of half-filled zero- 
energy Landau level, which is eightfold degenerate due 
to flavor (spin and valley (layer)) and orbital (n = 0, 1 
Landau index) degrees of freedom. We expect a gradual 
formation of a QHFM (quantum Hall ferromagnet) [l^ 
due to many-body correlations and the spontaneous fer- 
romagnetic ordering of the spin degree of freedom. Thus 
we will fix the valley and orbital degree of freedom in 
the following and discuss how from two (spin up and 
spin down) Landau levels we can have effectively a single 



filled Landau level and ferromagnetic ordering, d + id 
wave Cooper pairs, described in the long-distance with 
the following Cooper pair wave function 37 1 



(23) 



in the presence of the magnetic flux will be modified by 
fiux (vortex) attachment due to the particles of opposite 
spin as in 



/(n-f - r2i) 



Zi-\ — Z2l 
Zl-f - Z2l 



Y[{zit - Ziif 'W{z2i - Zjt) 



This will lead to the following many-body state 



(24) 



Det{ 



) n^^'t - - n^^'t - ^u)x2 (25) 



^,3 



where Det denotes the determinant of antisymmetrized 
product of Cooper pair wave functions, and X2 denotes 
the filled second Landau level wave function in the Jain 
notation. The identity used in Ea. (|25p was proved in 
Ref. [3^ The topological properties of the wave function 
in Eq. ipS]) (or the low energy properties of the system 
described with the wave function as discussed in Ref. [39l ) 
are equivalent to the Halperin (1,1,1) state or QHFM i.e. 
the following state 

J|(z,t - Zjt) Y[iZpl - Zql) Ylizi^ - Z„^l), (26) 
i<j p<q l,m 

for fixed valley and orbital index, and thus lead to the 
QHFM state with the effective filling factor i^off = 4. It 
was shown in Ref. that this state would lead to the 
gap with linear dependence on the (strong) magnetic field 
as observed in Ref. [l]. Thus we described a possible route 
from antiferromagnetic state with d + id superconducting 
correlations to the spin QHFM state consistent with the 
experiment. 

M.V.M. thanks M. Goerbig, D. Tanaskovic, and J. 
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supported by the Serbian Ministry of Education and Sci- 
ence under project No. ON171017. 
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Energy minimization 

To find whether d—id or d + id SC correlations coexist 
with antiferromagnetism in the presence of the magnetic 
field, which direction is defined as in Eq.([T8|), we should 
compare the two ground state energies, 



E' 



d—id 



^[-TO - (n + l)(n + 2)^] + 

n" 

Y^[-E-n{n-l)5], 



(27) 



and 



E 



d-\-id 



n=0 
n" 

[— m — n{n — 1)(5] , 



(28) 



where ^ = 5^, and the bounds for the summations are 
determined by the lower cut-off, —E^, i.e. we have, in 
the d — id case, 



(n'-f 2)(n' + l) 



Er — m 



and 



and 



and 



n"{n" - 1) 



Er-E 



{n' + 2){n' + 1) 



Ec-E 



n"{n" -l) = 



Er — m 



(29) 



(30) 



(31) 



(32) 



in the d + id case. After a few steps of simple algebra we 
get 



E' 



■d—id _ ^d-\-id 



2{m- E). 



(33) 



Because E > m, the energy minimization favors d—id 
SC correlations for the fixed direction of the magnetic 
field (Eq.dlHl)). 



